home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / docs / apps / database / postgres / postgre4.z / postgre4 / src / utils / adt / geo-ops.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-04-21  |  43.0 KB  |  2,153 lines
  1. /*
  2.  * geo-ops.c --
  3.  *    2D geometric operations
  4.  */
  5.  
  6. #ifndef lint
  7. static char rcs_id[] = 
  8. "$Header: /private/postgres/src/utils/adt/RCS/geo-ops.c,v 1.27 1992/08/25 17:50:23 mer Exp $";
  9. #endif
  10.  
  11. #include <stdio.h>
  12. #include <strings.h>
  13. #include <math.h>
  14.  
  15. #include "utils/geo-decls.h"
  16. #include "utils/log.h"
  17.  
  18. #define LDELIM        '('
  19. #define RDELIM        ')'
  20. #define    DELIM        ','
  21. #define BOXNARGS    4
  22. #define    LSEGNARGS    4
  23. #define    POINTNARGS    2
  24.  
  25. #ifdef sequent
  26. #define HUGE_VAL    1.8e+308
  27. #endif
  28.  
  29. /***********************************************************************
  30.  **
  31.  **     Routines for two-dimensional boxes.
  32.  **
  33.  ***********************************************************************/
  34.  
  35. /*----------------------------------------------------------
  36.  * Formatting and conversion routines.
  37.  *---------------------------------------------------------*/
  38.  
  39. /*    box_in    -    convert a string to internal form.
  40.  *
  41.  *    str:    input string "(f8, f8, f8, f8)"
  42.  */
  43. BOX *
  44. box_in(str)
  45.     char    *str;
  46. {
  47.     double    tmp;
  48.     char    *p, *coord[BOXNARGS];
  49.     int    i;
  50.     BOX    *result;
  51.     double    atof();
  52.  
  53.     if (str == NULL)
  54.         return(NULL);
  55.     for (i = 0, p = str; *p && i < BOXNARGS && *p != RDELIM; p++)
  56.         if (*p == DELIM || (*p == LDELIM && !i))
  57.             coord[i++] = p + 1;
  58.     if (i < BOXNARGS - 1)
  59.         return(NULL);
  60.     result = PALLOCTYPE(BOX);
  61.     result->xh = atof(coord[0]);
  62.     result->yh = atof(coord[1]);
  63.     result->xl = atof(coord[2]);
  64.     result->yl = atof(coord[3]);
  65.     if (result->xh < result->xl) {
  66.         tmp = result->xh;
  67.         result->xh = result->xl;
  68.         result->xl = tmp;
  69.     }
  70.     if (result->yh < result->yl) {
  71.         tmp = result->yh;
  72.         result->yh = result->yl;
  73.         result->yl = tmp;
  74.     }
  75.  
  76.     return(result);
  77. }
  78.  
  79. /*    box_out    -    convert a box to external form.
  80.  */
  81. char *
  82. box_out(box)
  83.     BOX    *box;
  84. {
  85.     char    *result;
  86.  
  87.     if (box == NULL)
  88.         return(NULL);
  89.     result = (char *)PALLOC(80);
  90.     (void) sprintf(result, "(%g,%g,%g,%g)",
  91.                box->xh, box->yh, box->xl, box->yl);
  92.  
  93.     return(result);
  94. }
  95.  
  96.  
  97. /*    box_construct    -    fill in a new box.
  98.  */
  99. BOX *
  100. box_construct(x1, x2, y1, y2)
  101.     double    x1, x2, y1, y2;
  102. {
  103.     BOX    *result;
  104.  
  105.     result = PALLOCTYPE(BOX);
  106.     return( box_fill(result, x1, x2, y1, y2) );
  107. }
  108.  
  109.  
  110. /*    box_fill    -    fill in a static box
  111.  */
  112. BOX *
  113. box_fill(result, x1, x2, y1, y2)
  114.     BOX    *result;
  115.     double    x1, x2, y1, y2;
  116. {
  117.     double    tmp;
  118.  
  119.     result->xh = x1;
  120.     result->xl = x2;
  121.     result->yh = y1;
  122.     result->yl = y2;
  123.     if (result->xh < result->xl) {
  124.         tmp = result->xh;
  125.         result->xh = result->xl;
  126.         result->xl = tmp;
  127.     }
  128.     if (result->yh < result->yl) {
  129.         tmp = result->yh;
  130.         result->yh = result->yl;
  131.         result->yl = tmp;
  132.     }
  133.  
  134.     return(result);
  135. }
  136.  
  137.  
  138. /*    box_copy    -    copy a box
  139.  */
  140. BOX *
  141. box_copy(box)
  142.     BOX    *box;
  143. {
  144.     BOX    *result;
  145.  
  146.     result = PALLOCTYPE(BOX);
  147.     bcopy((char *) box, (char *) result, sizeof(BOX));
  148.  
  149.     return(result);
  150. }
  151.  
  152.  
  153. /*----------------------------------------------------------
  154.  *  Relational operators for BOXes.
  155.  *    <, >, <=, >=, and == are based on box area.
  156.  *---------------------------------------------------------*/
  157.  
  158. /*    box_same    -    are two boxes identical?
  159.  */
  160. long
  161. box_same(box1, box2)
  162.     BOX    *box1, *box2;
  163. {
  164.     return((box1->xh == box2->xh && box1->xl == box2->xl) &&
  165.            (box1->yh == box2->yh && box1->yl == box2->yl));
  166. }
  167.  
  168. /*    box_overlap    -    does box1 overlap box2?
  169.  */
  170. long
  171. box_overlap(box1, box2)
  172.     BOX    *box1, *box2;
  173. {
  174.     return(((box1->xh >= box2->xh && box1->xl <= box2->xh) ||
  175.         (box2->xh >= box1->xh && box2->xl <= box1->xh)) &&
  176.            ((box1->yh >= box2->yh && box1->yl <= box2->yh) ||
  177.         (box2->yh >= box1->yh && box2->yl <= box1->yh)) );
  178. }
  179.  
  180. /*    box_overleft    -    is the right edge of box1 to the left of
  181.  *                the right edge of box2?
  182.  *
  183.  *    This is "less than or equal" for the end of a time range,
  184.  *    when time ranges are stored as rectangles.
  185.  */
  186. long
  187. box_overleft(box1, box2)
  188.     BOX    *box1, *box2;
  189. {
  190.     return(box1->xh <= box2->xh);
  191. }
  192.  
  193. /*    box_left    -    is box1 strictly left of box2?
  194.  */
  195. long
  196. box_left(box1, box2)
  197.     BOX    *box1, *box2;
  198. {
  199.     return(box1->xh < box2->xl);
  200. }
  201.  
  202. /*    box_right    -    is box1 strictly right of box2?
  203.  */
  204. long
  205. box_right(box1, box2)
  206.     BOX    *box1, *box2;
  207. {
  208.     return(box1->xl > box2->xh);
  209. }
  210.  
  211. /*    box_overright    -    is the left edge of box1 to the right of
  212.  *                the left edge of box2?
  213.  *
  214.  *    This is "greater than or equal" for time ranges, when time ranges
  215.  *    are stored as rectangles.
  216.  */
  217. long
  218. box_overright(box1, box2)
  219.     BOX    *box1, *box2;
  220. {
  221.     return(box1->xl >= box2->xl);
  222. }
  223.  
  224. /*    box_contained    -    is box1 contained by box2?
  225.  */
  226. long
  227. box_contained(box1, box2)
  228.     BOX    *box1, *box2;
  229. {
  230.     return((box1->xh <= box2->xh && box1->xl >= box2->xl &&
  231.         box1->yh <= box2->yh && box1->yl >= box2->yl));
  232. }
  233.  
  234. /*    box_contain    -    does box1 contain box2?
  235.  */
  236. long
  237. box_contain(box1, box2)
  238.     BOX    *box1, *box2;
  239. {
  240.     return((box1->xh >= box2->xh && box1->xl <= box2->xl &&
  241.         box1->yh >= box2->yh && box1->yl <= box2->yl));
  242. }
  243.  
  244.  
  245. /*    box_positionop    -
  246.  *        is box1 entirely {above, below } box2?
  247.  */
  248. long
  249. box_below(box1, box2)
  250.     BOX    *box1, *box2;
  251. { return( box1->yh <= box2->yl ); }
  252.  
  253. long
  254. box_above(box1, box2)
  255.     BOX    *box1, *box2;
  256. { return( box1->yl >= box2->yh ); }
  257.  
  258.  
  259. /*    box_relop    -    is area(box1) relop area(box2), within
  260.  *                  our accuracy constraint?
  261.  */
  262. long
  263. box_lt(box1, box2)
  264.     BOX    *box1, *box2;
  265. {
  266.     return( FPlt(box_ar(box1), box_ar(box2)) );
  267. }
  268.  
  269. long
  270. box_gt(box1, box2)
  271.     BOX    *box1, *box2;
  272. {
  273.     return( FPgt(box_ar(box1), box_ar(box2)) );
  274. }
  275.  
  276. long
  277. box_eq(box1, box2)
  278.     BOX    *box1, *box2;
  279. {
  280.     return( FPeq(box_ar(box1), box_ar(box2)) );
  281. }
  282.  
  283. long
  284. box_le(box1, box2)
  285.     BOX    *box1, *box2;
  286. {
  287.     return( FPle(box_ar(box1), box_ar(box2)) );
  288. }
  289.  
  290. long
  291. box_ge(box1, box2)
  292.     BOX    *box1, *box2;
  293. {
  294.     return( FPge(box_ar(box1), box_ar(box2)) );
  295. }
  296.  
  297.  
  298. /*----------------------------------------------------------
  299.  *  "Arithmetic" operators on boxes.
  300.  *    box_foo    returns foo as an object (pointer) that
  301.  can be passed between languages.
  302.  *    box_xx    is an internal routine which returns the
  303.  *        actual value (and cannot be handed back to
  304.  *        LISP).
  305.  *---------------------------------------------------------*/
  306.  
  307. /*    box_area    -    returns the area of the box.
  308.  */
  309. double *
  310. box_area(box)
  311.     BOX    *box;
  312. {
  313.     double *result;
  314.  
  315.     result = PALLOCTYPE(double);
  316.     *result = box_ln(box) * box_ht(box);
  317.  
  318.     return(result);
  319. }
  320.  
  321.  
  322. /*    box_length    -    returns the length of the box 
  323.  *                  (horizontal magnitude).
  324.  */
  325. double *
  326. box_length(box)
  327.     BOX    *box;
  328. {
  329.     double    *result;
  330.  
  331.     result = PALLOCTYPE(double);
  332.     *result = box->xh - box->xl;
  333.  
  334.     return(result);
  335. }
  336.  
  337.  
  338. /*    box_height    -    returns the height of the box 
  339.  *                  (vertical magnitude).
  340.  */
  341. double *
  342. box_height(box)
  343.     BOX    *box;
  344. {
  345.     double    *result;
  346.  
  347.     result = PALLOCTYPE(double);
  348.     *result = box->yh - box->yl;
  349.  
  350.     return(result);
  351. }
  352.  
  353.  
  354. /*    box_distance    -    returns the distance between the
  355.  *                  center points of two boxes.
  356.  */
  357. double *
  358. box_distance(box1, box2)
  359.     BOX    *box1, *box2;
  360. {
  361.     double    *result;
  362.     POINT    *box_center(), *a, *b;
  363.  
  364.     result = PALLOCTYPE(double);
  365.     a = box_center(box1);
  366.     b = box_center(box2);
  367.     *result = HYPOT(a->x - b->x, a->y - b->y);
  368.  
  369.     PFREE(a);
  370.     PFREE(b);
  371.     return(result);
  372. }
  373.  
  374.  
  375. /*    box_center    -    returns the center point of the box.
  376.  */
  377. POINT *
  378. box_center(box)
  379.     BOX    *box;
  380. {
  381.     POINT    *result;
  382.  
  383.     result = PALLOCTYPE(POINT);
  384.     result->x = (box->xh + box->xl) / 2.0;
  385.     result->y = (box->yh + box->yl) / 2.0;
  386.  
  387.     return(result);
  388. }
  389.  
  390.  
  391. /*    box_ar    -    returns the area of the box.
  392.  */
  393. double
  394. box_ar(box)
  395.     BOX    *box;
  396. {
  397.     return( box_ln(box) * box_ht(box) );
  398. }
  399.  
  400.  
  401. /*    box_ln    -    returns the length of the box 
  402.  *                  (horizontal magnitude).
  403.  */
  404. double
  405. box_ln(box)
  406.     BOX    *box;
  407. {
  408.     return( box->xh - box->xl );
  409. }
  410.  
  411.  
  412. /*    box_ht    -    returns the height of the box 
  413.  *                  (vertical magnitude).
  414.  */
  415. double
  416. box_ht(box)
  417.     BOX    *box;
  418. {
  419.     return( box->yh - box->yl );
  420. }
  421.  
  422.  
  423. /*    box_dt    -    returns the distance between the
  424.  *              center points of two boxes.
  425.  */
  426. double
  427. box_dt(box1, box2)
  428.     BOX    *box1, *box2;
  429. {
  430.     double    result;
  431.     POINT    *box_center(),
  432.     *a, *b;
  433.  
  434.     a = box_center(box1);
  435.     b = box_center(box2);
  436.     result = HYPOT(a->x - b->x, a->y - b->y);
  437.  
  438.     PFREE(a);
  439.     PFREE(b);
  440.     return(result);
  441. }
  442.  
  443. /*----------------------------------------------------------
  444.  *  Funky operations.
  445.  *---------------------------------------------------------*/
  446.  
  447. /*    box_intersect    -
  448.  *        returns the overlapping portion of two boxes,
  449.  *          or NULL if they do not intersect.
  450.  */
  451. BOX *
  452. box_intersect(box1, box2)
  453.     BOX    *box1, *box2;
  454. {
  455.     BOX    *result;
  456.     long    box_overlap();
  457.  
  458.     if (! box_overlap(box1,box2))
  459.         return(NULL);
  460.     result = PALLOCTYPE(BOX);
  461.     result->xh = MIN(box1->xh, box2->xh);
  462.     result->xl = MAX(box1->xl, box2->xl);
  463.     result->yh = MIN(box1->yh, box2->yh);
  464.     result->yl = MAX(box1->yl, box2->yl);
  465.  
  466.     return(result);
  467. }
  468.  
  469.  
  470. /*    box_diagonal    -    
  471.  *        returns a line segment which happens to be the
  472.  *          positive-slope diagonal of "box".
  473.  *        provided, of course, we have LSEGs.
  474.  */
  475. LSEG *
  476. box_diagonal(box)
  477.     BOX    *box;
  478. {
  479.     POINT    p1, p2;
  480.  
  481.     p1.x = box->xh;
  482.     p1.y = box->yh;
  483.     p2.x = box->xl;
  484.     p2.y = box->yl;
  485.     return( lseg_construct( &p1, &p2 ) );
  486.  
  487. }
  488.  
  489. /***********************************************************************
  490.  **
  491.  **     Routines for 2D lines.
  492.  **        Lines are not intended to be used as ADTs per se,
  493.  **        but their ops are useful tools for other ADT ops.  Thus,
  494.  **        there are few relops.
  495.  **
  496.  ***********************************************************************/
  497.  
  498. /*----------------------------------------------------------
  499.  *  Conversion routines from one line formula to internal.
  500.  *    Internal form:    Ax+By+C=0
  501.  *---------------------------------------------------------*/
  502.  
  503. LINE *                /* point-slope */
  504. line_construct_pm(pt, m)
  505.     POINT    *pt;
  506.     double    m;
  507. {
  508.     LINE    *result;
  509.  
  510.     result = PALLOCTYPE(LINE);
  511.     /* use "mx - y + yinter = 0" */
  512.     result->A = m;
  513.     result->B = -1.0;
  514.     result->C = pt->y - m * pt->x;
  515.     return(result);
  516. }
  517.  
  518.  
  519. LINE *                /* two points */
  520. line_construct_pp(pt1, pt2)
  521.     POINT    *pt1, *pt2;
  522. {
  523.     LINE    *result;
  524.  
  525.     result = PALLOCTYPE(LINE);
  526.     if (FPeq(pt1->x, pt2->x)) {        /* vertical */
  527.         /* use "x = C" */
  528.         result->m = 0.0;
  529.         result->A = -1.0;
  530.         result->B = 0.0;
  531.         result->C = pt1->x;
  532.     } else {
  533.         /* use "mx - y + yinter = 0" */
  534.         result->m = (pt1->y - pt2->y) / (pt1->x - pt2->x);
  535.         result->A = result->m;
  536.         result->B = -1.0;
  537.         result->C = pt1->y - result->m * pt1->x;
  538.     }
  539.     return(result);
  540. }
  541.  
  542.  
  543. /*----------------------------------------------------------
  544.  *  Relative position routines.
  545.  *---------------------------------------------------------*/
  546.  
  547. long
  548. line_intersect(l1, l2)
  549.     LINE    *l1, *l2;
  550. {
  551.     return( ! line_parallel(l1, l2) );
  552. }
  553.  
  554. long
  555. line_parallel(l1, l2)
  556.     LINE    *l1, *l2;
  557. {
  558.     return( FPeq(l1->m, l2->m) );
  559. }
  560.  
  561. long
  562. line_perp(l1, l2)
  563.     LINE    *l1, *l2;
  564. {
  565.     if (l1->m)
  566.         return( FPeq(l2->m / l1->m, -1.0) );
  567.     else if (l2->m)
  568.         return( FPeq(l1->m / l2->m, -1.0) );
  569.     return(1);    /* both 0.0 */
  570. }
  571.  
  572. long
  573. line_vertical(line)
  574.     LINE    *line;
  575. {
  576.     return( FPeq(line->A, -1.0) && FPzero(line->B) );
  577. }
  578.  
  579. long
  580. line_horizontal(line)
  581.     LINE    *line;
  582. {
  583.     return( FPzero(line->m) );
  584. }
  585.  
  586.  
  587. long
  588. line_eq(l1, l2)
  589.     LINE    *l1, *l2;
  590. {
  591.     double    k;
  592.  
  593.     if (! FPzero(l2->A))
  594.         k = l1->A / l2->A;
  595.     else if (! FPzero(l2->B))
  596.         k = l1->B / l2->B;
  597.     else if (! FPzero(l2->C))
  598.         k = l1->C / l2->C;
  599.     else
  600.         k = 1.0;
  601.     return( FPeq(l1->A, k * l2->A) &&
  602.            FPeq(l1->B, k * l2->B) &&
  603.            FPeq(l1->C, k * l2->C) );
  604. }
  605.  
  606.  
  607. /*----------------------------------------------------------
  608.  *  Line arithmetic routines.
  609.  *---------------------------------------------------------*/
  610.  
  611. double *         /* distance between l1, l2 */
  612. line_distance(l1, l2)
  613.     LINE    *l1, *l2;
  614. {
  615.     double    *result;
  616.     POINT    *tmp;
  617.  
  618.     result = PALLOCTYPE(double);
  619.     if (line_intersect(l1, l2)) {
  620.         *result = 0.0;
  621.         return(result);
  622.     }
  623.     if (line_vertical(l1))
  624.         *result = fabs(l1->C - l2->C);
  625.     else {
  626.         tmp = point_construct(0.0, l1->C);
  627.         result = dist_pl(tmp, l2);
  628.         PFREE(tmp);
  629.     }
  630.     return(result);
  631. }
  632.  
  633. POINT *            /* point where l1, l2 intersect (if any) */
  634.     line_interpt(l1, l2)
  635.     LINE    *l1, *l2;
  636. {
  637.     POINT    *result;
  638.     double    x;
  639.  
  640.     if (line_parallel(l1, l2))
  641.         return(NULL);
  642.     if (line_vertical(l1))
  643.         result = point_construct(l2->m * l1->C + l2->C, l1->C);
  644.     else if (line_vertical(l2))
  645.         result = point_construct(l1->m * l2->C + l1->C, l2->C);
  646.     else {
  647.         x = (l1->C - l2->C) / (l2->A - l1->A);
  648.         result = point_construct(x, l1->m * x + l1->C);
  649.     }
  650.     return(result);
  651. }
  652.  
  653. /***********************************************************************
  654.  **
  655.  **     Routines for 2D paths (sequences of line segments, also
  656.  **        called `polylines').
  657.  **
  658.  **        This is not a general package for geometric paths, 
  659.  **        which of course include polygons; the emphasis here
  660.  **        is on (for example) usefulness in wire layout.
  661.  **
  662.  ***********************************************************************/
  663.  
  664. #define    PATHALLOCSIZE(N) \
  665. (long) ((unsigned) (sizeof(PATH) + \
  666.                 (((N)-1) > 0 ? ((N)-1) : 0) \
  667.                 * sizeof(POINT)))
  668.  
  669. /*----------------------------------------------------------
  670.  *  String to path / path to string conversion.
  671.  *    External format: 
  672.  *        "(closed, npts, xcoord, ycoord,... )"
  673.  *---------------------------------------------------------*/
  674.  
  675. PATH *
  676. path_in(str)
  677.     char    *str;
  678. {
  679.     double    atof(), coord;
  680.     long    atol(), field[2];
  681.     char    *s;
  682.     int    ct, i;
  683.     PATH    *result;
  684.     long    pathsize;
  685.  
  686.     if (str == NULL)
  687.         elog(WARN, "Bad (null) path string representation");
  688.  
  689.     /* read the path header information */
  690.     for (i = 0, s = str; *s && i < 2 && *s != RDELIM; ++s)
  691.         if (*s == DELIM || (*s == LDELIM && !i))
  692.             field[i++] = atol(s + 1);
  693.     if (i < 1)
  694.         elog(WARN, "Bad path string representation: %s", str);
  695.     pathsize = PATHALLOCSIZE(field[1]);
  696.     result = (PATH *)palloc(pathsize);
  697.     result->length = pathsize;
  698.     result->closed = field[0];
  699.     result->npts =  field[1];
  700.  
  701.     /* read the path points */
  702.  
  703.     ct = result->npts * 2;    /* two coords for every point */
  704.     for (i = 0;
  705.          *s && i < ct && *s != RDELIM; 
  706.          ++s) {
  707.         if (*s == ',') {
  708.             coord = atof(s + 1);
  709.             if (i % 2)
  710.                 (result->p[i/2]).y = coord;
  711.             else
  712.                 (result->p[i/2]).x = coord;
  713.             ++i;
  714.         }
  715.     }
  716.     if (i % 2 || i < --ct) {
  717.         PFREE(result);
  718.         elog(WARN, "Bad path string representation: %s", str);
  719.     } 
  720.  
  721.     return(result);
  722. }
  723.  
  724.  
  725. char *
  726. path_out(path)
  727.     PATH    *path;
  728. {
  729.     char         buf[BUFSIZ + 20000], *result, *s;
  730.     int        i;
  731.     char    tmp[64];
  732.  
  733.     if (path == NULL)
  734.         return(NULL);
  735.     (void) sprintf(buf,"%c%d,%d\0", LDELIM, 
  736.                path->closed, path->npts);
  737.     s = buf + strlen(buf);
  738.     for (i = 0; i < path->npts; ++i) {
  739.         (void) sprintf(tmp, ",%g,%g", 
  740.                    path->p[i].x, path->p[i].y);
  741.         (void) strcpy(s, tmp);
  742.         s += strlen(tmp);
  743.     }
  744.     *s++ = RDELIM;
  745.     *s = '\0';
  746.     result = (char *)PALLOC(strlen(buf) + 1);
  747.     (void) strcpy(result, buf);
  748.  
  749.     return(result);
  750. }
  751.  
  752.  
  753. /*----------------------------------------------------------
  754.  *  Relational operators.
  755.  *    These are based on the path cardinality, 
  756.  *    as stupid as that sounds.
  757.  *
  758.  *    Better relops and access methods coming soon.
  759.  *---------------------------------------------------------*/
  760.  
  761. long
  762. path_n_lt(p1, p2)
  763.     PATH     *p1, *p2;
  764. {
  765.     return( (p1->npts < p2->npts ) );
  766. }
  767.  
  768. long
  769. path_n_gt(p1, p2)
  770.     PATH     *p1, *p2;
  771. {
  772.     return( (p1->npts > p2->npts ) );
  773. }
  774.  
  775. long
  776. path_n_eq(p1, p2)
  777.     PATH     *p1, *p2;
  778. {
  779.     return( (p1->npts == p2->npts) );
  780. }
  781.  
  782. long
  783. path_n_le(p1, p2)
  784.     PATH     *p1, *p2;
  785. {
  786.     return( (p1->npts <= p2->npts ) );
  787. }
  788.  
  789. long
  790. path_n_ge(p1, p2)
  791.     PATH     *p1, *p2;
  792. {
  793.     return( (p1->npts >= p2->npts ) );
  794. }
  795.  
  796. /* path_inter -
  797.  *    Does p1 intersect p2 at any point?
  798.  *    Use bounding boxes for a quick (O(n)) check, then do a 
  799.  *    O(n^2) iterative edge check.
  800.  */
  801. long
  802. path_inter(p1, p2)
  803.      PATH    *p1, *p2;
  804. {
  805.     BOX    b1, b2;
  806.     int    i, j;
  807.     LSEG seg1, seg2;
  808.     
  809.     b1.xh = b1.yh = b2.xh = b2.yh = HUGE;
  810.     b1.xl = b1.yl = b2.xl = b2.yl = -HUGE;
  811.     for (i = 0; i < p1->npts; ++i) {
  812.     b1.xh = MAX(p1->p[i].x, b1.xh);
  813.     b1.yh = MAX(p1->p[i].y, b1.yh);
  814.     b1.xl = MIN(p1->p[i].x, b1.xl);
  815.     b1.yl = MIN(p1->p[i].y, b1.yl);
  816.     }
  817.     for (i = 0; i < p2->npts; ++i) {
  818.     b2.xh = MAX(p2->p[i].x, b2.xh);
  819.     b2.yh = MAX(p2->p[i].y, b2.yh);
  820.     b2.xl = MIN(p2->p[i].x, b2.xl);
  821.     b2.yl = MIN(p2->p[i].y, b2.yl);
  822.     }
  823.     if (! box_overlap(&b1, &b2))
  824.       return(0);
  825.     
  826.     /*  pairwise check lseg intersections */
  827.     for (i = 0; i < p1->npts - 1; i++)
  828.       for (j = 0; j < p2->npts - 1; j++)
  829.        {
  830.        statlseg_construct(&seg1, &p1->p[i], &p1->p[i+1]);
  831.        statlseg_construct(&seg2, &p2->p[j], &p2->p[j+1]);
  832.        if (lseg_intersect(&seg1, &seg2))
  833.          return(1);
  834.        }
  835.     
  836.     /* if we dropped through, no two segs intersected */
  837.     return(0);
  838. }
  839.  
  840. /* this essentially does a cartesian product of the lsegs in the
  841.    two paths, and finds the min distance between any two lsegs */
  842. double *path_distance(p1, p2)
  843.      PATH *p1;
  844.      PATH *p2;
  845. {
  846.     double *min, *tmp;
  847.     int i,j;
  848.     LSEG seg1, seg2;
  849.  
  850.     statlseg_construct(&seg1, &p1->p[0], &p1->p[1]);
  851.     statlseg_construct(&seg2, &p2->p[0], &p2->p[1]);
  852.     min = lseg_distance(&seg1, &seg2);
  853.  
  854.     for (i = 0; i < p1->npts - 1; i++)
  855.       for (j = 0; j < p2->npts - 1; j++)
  856.        {
  857.        statlseg_construct(&seg1, &p1->p[i], &p1->p[i+1]);
  858.        statlseg_construct(&seg2, &p2->p[j], &p2->p[j+1]);
  859.  
  860.        if (*min < *(tmp = lseg_distance(&seg1, &seg2)))
  861.          *min = *tmp;
  862.        PFREE(tmp);
  863.        }
  864.  
  865.     return(min);
  866. }
  867.  
  868.  
  869. /*----------------------------------------------------------
  870.  *  "Arithmetic" operations.
  871.  *---------------------------------------------------------*/
  872.  
  873. double *
  874. path_length(path)
  875.     PATH    *path;
  876. {
  877.     double    *result;
  878.     int    ct, i;
  879.  
  880.     result = PALLOCTYPE(double);
  881.     ct = path->npts - 1;
  882.     for (i = 0; i < ct; ++i)
  883.         *result += point_dt(&path->p[i], &path->p[i+1]);
  884.  
  885.     return(result);
  886. }
  887.  
  888.  
  889.  
  890. double
  891. path_ln(path)
  892.     PATH    *path;
  893. {
  894.     double    result;
  895.     int    ct, i;
  896.  
  897.     ct = path->npts - 1;
  898.     for (result = i = 0; i < ct; ++i)
  899.         result += point_dt(&path->p[i], &path->p[i+1]);
  900.  
  901.     return(result);
  902. }
  903. /***********************************************************************
  904.  **
  905.  **     Routines for 2D points.
  906.  **
  907.  ***********************************************************************/
  908.  
  909. /*----------------------------------------------------------
  910.  *  String to point, point to string conversion.
  911.  *    External form:    "(x, y)"
  912.  *---------------------------------------------------------*/
  913.  
  914. POINT *
  915. point_in(str)
  916.     char    *str;
  917. {
  918.     double    atof();
  919.     char    *coord[POINTNARGS], *p;
  920.     int    i;
  921.     POINT    *result;
  922.  
  923.     if (str == NULL)
  924.         elog(WARN, "Bad (Null) point external representation");
  925.     for (i = 0, p = str; *p && i < POINTNARGS && *p != RDELIM; p++)
  926.         if (*p == DELIM || (*p == LDELIM && !i))
  927.             coord[i++] = p + 1;
  928.     if (i < POINTNARGS - 1)
  929.         elog(WARN, "Bad point external representation '%s'",str);
  930.     result = PALLOCTYPE(POINT);
  931.     result->x = atof(coord[0]);
  932.     result->y = atof(coord[1]);
  933.     return(result);
  934. }
  935.  
  936.  
  937. char *
  938. point_out(pt)
  939.     POINT    *pt;
  940. {
  941.     char    *result;
  942.  
  943.     if (pt == NULL)
  944.         return(NULL);
  945.     result = (char *)PALLOC(40);
  946.     (void) sprintf(result, "(%g,%g)", pt->x, pt->y);
  947.     return(result);
  948. }
  949.  
  950.  
  951. POINT *
  952. point_construct(x, y)
  953.     double    x, y;
  954. {
  955.     POINT    *result;
  956.  
  957.     result = PALLOCTYPE(POINT);
  958.     result->x = x;
  959.     result->y = y;
  960.     return(result);
  961. }
  962.  
  963.  
  964. POINT *
  965. point_copy(pt)
  966.     POINT    *pt;
  967. {
  968.     POINT    *result;
  969.  
  970.     result = PALLOCTYPE(POINT);
  971.     result->x = pt->x;
  972.     result->y = pt->y;
  973.     return(result);
  974. }
  975.  
  976.  
  977. /*----------------------------------------------------------
  978.  *  Relational operators for POINTs.
  979.  *    Since we do have a sense of coordinates being
  980.  *    "equal" to a given accuracy (point_vert, point_horiz), 
  981.  *    the other ops must preserve that sense.  This means
  982.  *    that results may, strictly speaking, be a lie (unless
  983.  *    EPSILON = 0.0).
  984.  *---------------------------------------------------------*/
  985.  
  986. long
  987. point_left(pt1, pt2)
  988.     POINT    *pt1, *pt2;
  989. { return( FPlt(pt1->x, pt2->x) ); }
  990.  
  991. long
  992. point_right(pt1, pt2)
  993.     POINT    *pt1, *pt2;
  994. { return( FPgt(pt1->x, pt2->x) ); }
  995.  
  996. long
  997. point_above(pt1, pt2)
  998.     POINT    *pt1, *pt2;
  999. { return( FPgt(pt1->y, pt2->y) ); }
  1000.  
  1001. long
  1002. point_below(pt1, pt2)
  1003.     POINT    *pt1, *pt2;
  1004. { return( FPlt(pt1->y, pt2->y) ); }
  1005.  
  1006. long
  1007. point_vert(pt1, pt2)
  1008.     POINT    *pt1, *pt2;
  1009. { return( FPeq( pt1->x, pt2->x ) ); }
  1010.  
  1011. long
  1012. point_horiz(pt1, pt2)
  1013.     POINT    *pt1, *pt2;
  1014. { return( FPeq( pt1->y, pt2->y ) ); }
  1015.  
  1016. long
  1017. point_eq(pt1, pt2)
  1018.     POINT    *pt1, *pt2;
  1019. { return( point_horiz(pt1, pt2) && point_vert(pt1, pt2) ); }
  1020.  
  1021. /*----------------------------------------------------------
  1022.  *  "Arithmetic" operators on points.
  1023.  *---------------------------------------------------------*/
  1024.  
  1025. long
  1026. pointdist(p1, p2)
  1027.     POINT    *p1, *p2;
  1028. {
  1029.     long result;
  1030.  
  1031.     result = point_dt(p1, p2);
  1032.     return(result);
  1033. }
  1034.  
  1035. double *
  1036. point_distance(pt1, pt2)
  1037.     POINT    *pt1, *pt2;
  1038. {
  1039.     double    *result;
  1040.  
  1041.     result = PALLOCTYPE(double);
  1042.     *result = HYPOT( pt1->x - pt2->x, pt1->y - pt2->y );
  1043.     return(result);
  1044. }
  1045.  
  1046.  
  1047. double
  1048. point_dt(pt1, pt2)
  1049.     POINT    *pt1, *pt2;
  1050. {
  1051.     return( HYPOT( pt1->x - pt2->x, pt1->y - pt2->y ) );
  1052. }
  1053.  
  1054. double *
  1055. point_slope(pt1, pt2)
  1056.     POINT    *pt1, *pt2;
  1057. {
  1058.     double    *result;
  1059.  
  1060.     result = PALLOCTYPE(double);
  1061.     if (point_vert(pt1, pt2))
  1062.         *result = HUGE;
  1063.     else
  1064.         *result = (pt1->y - pt2->y) / (pt1->x - pt1->x);
  1065.     return(result);
  1066. }
  1067.  
  1068.  
  1069. double
  1070. point_sl(pt1, pt2)
  1071.     POINT    *pt1, *pt2;
  1072. {
  1073.     return(    point_vert(pt1, pt2)
  1074.            ? HUGE
  1075.            : (pt1->y - pt2->y) / (pt1->x - pt2->x) );
  1076. }
  1077.  
  1078. /***********************************************************************
  1079.  **
  1080.  **     Routines for 2D line segments.
  1081.  **
  1082.  ***********************************************************************/
  1083.  
  1084. /*----------------------------------------------------------
  1085.  *  String to lseg, lseg to string conversion.
  1086.  *    External form:    "(id, info, x1, y1, x2, y2)"
  1087.  *---------------------------------------------------------*/
  1088.  
  1089. LSEG *
  1090. lseg_in(str)
  1091.     char    *str;
  1092. {
  1093.     char    *coord[LSEGNARGS], *p;
  1094.     int    i;
  1095.     LSEG    *result;
  1096.     extern double    atof();
  1097.  
  1098.     if (str == NULL)
  1099.         return(NULL);
  1100.     for (i = 0, p = str; *p && i < LSEGNARGS && *p != RDELIM; p++)
  1101.         if (*p == DELIM || (*p == LDELIM && !i))
  1102.             coord[i++] = p + 1;
  1103.     if (i < LSEGNARGS - 1)
  1104.         return(NULL);
  1105.     result = PALLOCTYPE(LSEG);
  1106.     result->p[0].x = atof(coord[0]);
  1107.     result->p[0].y = atof(coord[1]);
  1108.     result->p[1].x = atof(coord[2]);
  1109.     result->p[1].y = atof(coord[3]);
  1110.     result->m = point_sl(&result->p[0], &result->p[1]);
  1111.  
  1112.     return(result);
  1113. }
  1114.  
  1115.  
  1116. char *
  1117. lseg_out(ls)
  1118.     LSEG    *ls;
  1119. {
  1120.     char    *result;
  1121.  
  1122.     if (ls == NULL)
  1123.         return(NULL);
  1124.     result = (char *)PALLOC(80);
  1125.     (void) sprintf(result, "(%g,%g,%g,%g)",
  1126.                ls->p[0].x, ls->p[0].y, ls->p[1].x, ls->p[1].y);
  1127.     return(result);
  1128. }
  1129.  
  1130.  
  1131. /* lseg_construct -
  1132.  *    form a LSEG from two POINTs.
  1133.  */
  1134. LSEG *
  1135. lseg_construct(pt1, pt2)
  1136.     POINT    *pt1, *pt2;
  1137. {
  1138.     LSEG    *result;
  1139.  
  1140.     result = PALLOCTYPE(LSEG);
  1141.     result->p[0].x = pt1->x;
  1142.     result->p[0].y = pt1->y;
  1143.     result->p[1].x = pt2->x;
  1144.     result->p[1].y = pt2->y;
  1145.     result->m = point_sl(pt1, pt2);
  1146.  
  1147.     return(result);
  1148. }
  1149.  
  1150. /* like lseg_construct, but assume space already allocated */
  1151. void statlseg_construct(lseg, pt1, pt2)
  1152.      LSEG *lseg;
  1153.      POINT *pt1;
  1154.      POINT *pt2;
  1155. {
  1156.     lseg->p[0].x = pt1->x;
  1157.     lseg->p[0].y = pt1->y;
  1158.     lseg->p[1].x = pt2->x;
  1159.     lseg->p[1].y = pt2->y;
  1160.     lseg->m = point_sl(pt1, pt2);
  1161. }
  1162.  
  1163. /*----------------------------------------------------------
  1164.  *  Relative position routines.
  1165.  *---------------------------------------------------------*/
  1166.  
  1167. /*
  1168. **  find intersection of the two lines, and see if it falls on 
  1169. **  both segments.
  1170. */
  1171. long
  1172. lseg_intersect(l1, l2)
  1173.     LSEG    *l1, *l2;
  1174. {
  1175.     LINE *ln;
  1176.     POINT *interpt;
  1177.     long retval;
  1178.  
  1179.     ln = line_construct_pp(&l2->p[0], &l2->p[1]);
  1180.     interpt = interpt_sl(l1, ln);
  1181.  
  1182.     if (interpt != NULL && on_ps(interpt, l2)) /* interpt on l1 and l2 */
  1183.       retval = 1;
  1184.     else retval = 0;
  1185.     if (interpt != NULL) PFREE(interpt);
  1186.     PFREE(ln);
  1187.     return(retval);
  1188. }
  1189.  
  1190. long
  1191. lseg_parallel(l1, l2)
  1192.     LSEG    *l1, *l2;
  1193. {
  1194.     return( FPeq(l1->m, l2->m) );
  1195. }
  1196.  
  1197. long
  1198. lseg_perp(l1, l2)
  1199.     LSEG    *l1, *l2;
  1200. {
  1201.     if (! FPzero(l1->m))
  1202.         return( FPeq(l2->m / l1->m, -1.0) );
  1203.     else if (! FPzero(l2->m))
  1204.         return( FPeq(l1->m / l2->m, -1.0) );
  1205.     return(0);    /* both 0.0 */
  1206. }
  1207.  
  1208. long
  1209. lseg_vertical(lseg)
  1210.     LSEG    *lseg;
  1211. {
  1212.     return( FPeq(lseg->p[0].x, lseg->p[1].x) );
  1213. }
  1214.  
  1215. long
  1216. lseg_horizontal(lseg)
  1217.     LSEG    *lseg;
  1218. {
  1219.     return( FPeq(lseg->p[0].y, lseg->p[1].y) );
  1220. }
  1221.  
  1222.  
  1223. long
  1224. lseg_eq(l1, l2)
  1225.     LSEG    *l1, *l2;
  1226. {
  1227.     return( FPeq(l1->p[0].x, l2->p[0].x) &&
  1228.            FPeq(l1->p[1].y, l2->p[1].y) &&
  1229.            FPeq(l1->p[0].x, l2->p[0].x) &&
  1230.            FPeq(l1->p[1].y, l2->p[1].y) );
  1231. }
  1232.  
  1233.  
  1234. /*----------------------------------------------------------
  1235.  *  Line arithmetic routines.
  1236.  *---------------------------------------------------------*/
  1237.  
  1238. /* lseg_distance -
  1239.  *    If two segments don't intersect, then the closest
  1240.  *    point will be from one of the endpoints to the other
  1241.  *    segment.
  1242.  */
  1243. double *
  1244. lseg_distance(l1, l2)
  1245.     LSEG    *l1, *l2;
  1246. {
  1247.     double    *d, *result;
  1248.  
  1249.     result = PALLOCTYPE(double);
  1250.     if (lseg_intersect(l1, l2)) {
  1251.         *result = 0.0;
  1252.         return(result);
  1253.     }
  1254.     *result = HUGE;
  1255.     d = dist_ps(&l1->p[0], l2);
  1256.     *result = MIN(*result, *d);
  1257.     PFREE(d);
  1258.     d = dist_ps(&l1->p[1], l2);
  1259.     *result = MIN(*result, *d);
  1260.     PFREE(d);
  1261.     d = dist_ps(&l2->p[0], l1);
  1262.     *result = MIN(*result, *d);
  1263.     PFREE(d);
  1264.     d = dist_ps(&l2->p[1], l1);
  1265.     *result = MIN(*result, *d);
  1266.     PFREE(d);
  1267.  
  1268.     return(result);
  1269. }
  1270.  
  1271.  
  1272. double
  1273. lseg_dt(l1, l2)            /* distance between l1, l2 */
  1274.     LSEG    *l1, *l2;
  1275. {
  1276.     double    *d, result;
  1277.  
  1278.     if (lseg_intersect(l1, l2))
  1279.         return(0.0);
  1280.     result = HUGE;
  1281.     d = dist_ps(&l1->p[0], l2);
  1282.     result = MIN(result, *d);
  1283.     PFREE(d);
  1284.     d = dist_ps(&l1->p[1], l2);
  1285.     result = MIN(result, *d);
  1286.     PFREE(d);
  1287.     d = dist_ps(&l2->p[0], l1);
  1288.     result = MIN(result, *d);
  1289.     PFREE(d);
  1290.     d = dist_ps(&l2->p[1], l1);
  1291.     result = MIN(result, *d);
  1292.     PFREE(d);
  1293.  
  1294.     return(result);
  1295. }
  1296.  
  1297.  
  1298. /* lseg_interpt -
  1299.  *    Find the intersection point of two segments (if any).
  1300.  *    Find the intersection of the appropriate lines; if the 
  1301.  *    point is not on a given segment, there is no valid segment
  1302.  *    intersection point at all.
  1303.  */
  1304. POINT *
  1305. lseg_interpt(l1, l2)
  1306.     LSEG    *l1, *l2;
  1307. {
  1308.     POINT    *result;
  1309.     LINE    *tmp1, *tmp2;
  1310.  
  1311.     tmp1 = line_construct_pp(&l1->p[0], &l1->p[1]);
  1312.     tmp2 = line_construct_pp(&l2->p[0], &l2->p[1]);
  1313.     if (result = line_interpt(tmp1, tmp2))
  1314.         if (! on_ps(result, l1)) {
  1315.             PFREE(result);
  1316.             result = NULL;
  1317.         }
  1318.  
  1319.     PFREE(tmp1);
  1320.     PFREE(tmp2);
  1321.     return(result);
  1322. }
  1323.  
  1324. /***********************************************************************
  1325.  **
  1326.  **     Routines for position comparisons of differently-typed
  1327.  **        2D objects.
  1328.  **
  1329.  ***********************************************************************/
  1330.  
  1331. #define    ABOVE    1
  1332. #define    BELOW    0
  1333. #define    UNDEF    -1
  1334.  
  1335.  
  1336. /*---------------------------------------------------------------------
  1337.  *    dist_
  1338.  *        Minimum distance from one object to another.
  1339.  *-------------------------------------------------------------------*/
  1340.  
  1341. double *
  1342. dist_pl(pt, line)
  1343.     POINT    *pt;
  1344.     LINE    *line;
  1345. {
  1346.     double    *result;
  1347.  
  1348.     result = PALLOCTYPE(double);
  1349.     *result = (line->A * pt->x + line->B * pt->y + line->C) /
  1350.         HYPOT(line->A, line->B);
  1351.  
  1352.     return(result);
  1353. }
  1354.  
  1355. double *
  1356. dist_ps(pt, lseg)
  1357.      POINT *pt;
  1358.      LSEG *lseg;
  1359. {
  1360.     double m;                       /* slope of perp. */
  1361.     LINE *ln;
  1362.     double *result, *tmpdist;
  1363.     POINT *ip;
  1364.  
  1365.     /* construct a line that's perpendicular.  See if the intersection of
  1366.        the two lines is on the line segment. */
  1367.     if (lseg->p[1].x == lseg->p[0].x)
  1368.       m = 0;
  1369.     else if (lseg->p[1].y == lseg->p[0].y) /* slope is infinite */
  1370.       m = HUGE;
  1371.     else m = (-1) * (lseg->p[1].y - lseg->p[0].y) / 
  1372.                     (lseg->p[1].x - lseg->p[0].x);
  1373.     ln = line_construct_pm(pt, m);
  1374.  
  1375.     if ((ip = interpt_sl(lseg, ln)) != NULL)
  1376.       result = point_distance(pt, ip);
  1377.     else  /* intersection is not on line segment, so distance is min
  1378.          of distance from point to an endpoint */
  1379.      {
  1380.      result = point_distance(pt, &lseg->p[0]);
  1381.      tmpdist = point_distance(pt, &lseg->p[1]);
  1382.      if (*tmpdist < *result) *result = *tmpdist;
  1383.      PFREE (tmpdist);
  1384.      }
  1385.  
  1386.     if (ip != NULL) PFREE(ip);
  1387.     PFREE(ln);
  1388.     return (result);
  1389. }
  1390.  
  1391.  
  1392. /*
  1393. ** Distance from a point to a path 
  1394. */
  1395. double *dist_ppth(pt, path)
  1396.      POINT *pt;
  1397.      PATH *path;
  1398. {
  1399.     double *result;
  1400.     double *tmp;
  1401.     int i;
  1402.     LSEG lseg;
  1403.  
  1404.     if (path->npts = 0) 
  1405.      {
  1406.      result = PALLOCTYPE(double);
  1407.      *result = Abs((double) HUGE_VAL);
  1408.      goto exit;
  1409.      }
  1410.  
  1411.     if (path->npts = 1) 
  1412.       {
  1413.       result = point_distance(pt, &path->p[0]);
  1414.       goto exit;
  1415.       }
  1416.  
  1417.     /* else */
  1418.     for (i = 0; i < path->npts; i++)
  1419.      {
  1420.      statlseg_construct(&lseg, &path->p[i], &path->p[i+1]);
  1421.      if (i = 0) result = tmp = dist_ps(pt, &lseg);
  1422.      if (*tmp < *result) 
  1423.        *result = *tmp;
  1424.      PFREE(tmp);
  1425.      }
  1426.  
  1427.   exit:
  1428.     return(result);
  1429. }
  1430.  
  1431. double *
  1432. dist_pb(pt, box)
  1433.     POINT    *pt;
  1434.     BOX    *box;
  1435. {
  1436.     POINT    *tmp;
  1437.     double    *result;
  1438.  
  1439.     tmp = close_pb(pt, box);
  1440.     result = point_distance(tmp, pt);
  1441.  
  1442.     PFREE(tmp);
  1443.     return(result);
  1444. }
  1445.  
  1446.  
  1447. double *
  1448. dist_sl(lseg, line)
  1449.     LSEG    *lseg;
  1450.     LINE    *line;
  1451. {
  1452.     double    *result;
  1453.  
  1454.     if (inter_sl(lseg, line)) {
  1455.         result = PALLOCTYPE(double);
  1456.         *result = 0.0;
  1457.     } else    /* parallel */
  1458.         result = dist_pl(&lseg->p[0], line);
  1459.  
  1460.     return(result);
  1461. }
  1462.  
  1463.  
  1464. double *
  1465. dist_sb(lseg, box)
  1466.     LSEG    *lseg;
  1467.     BOX    *box;
  1468. {
  1469.     POINT    *tmp;
  1470.     double    *result;
  1471.  
  1472.     tmp = close_sb(lseg, box);
  1473.     if (tmp == NULL) {
  1474.         result = PALLOCTYPE(double);
  1475.         *result = 0.0;
  1476.     } else {
  1477.         result = dist_pb(tmp, box);
  1478.         PFREE(tmp);
  1479.     }
  1480.  
  1481.     return(result);
  1482. }
  1483.  
  1484.  
  1485. double *
  1486. dist_lb(line, box)
  1487.     LINE    *line;
  1488.     BOX    *box;
  1489. {
  1490.     POINT    *tmp;
  1491.     double    *result;
  1492.  
  1493.     tmp = close_lb(line, box);
  1494.     if (tmp == NULL) {
  1495.         result = PALLOCTYPE(double);
  1496.         *result = 0.0;
  1497.     } else {
  1498.         result = dist_pb(tmp, box);
  1499.         PFREE(tmp);
  1500.     }
  1501.  
  1502.     return(result);
  1503. }
  1504.  
  1505.  
  1506. /*---------------------------------------------------------------------
  1507.  *    interpt_
  1508.  *        Intersection point of objects.
  1509.  *        We choose to ignore the "point" of intersection between 
  1510.  *          lines and boxes, since there are typically two.
  1511.  *-------------------------------------------------------------------*/
  1512.  
  1513. POINT *
  1514. interpt_sl(lseg, line)
  1515.     LSEG    *lseg;
  1516.     LINE    *line;
  1517. {
  1518.     LINE    *tmp;
  1519.     POINT    *p;
  1520.  
  1521.     tmp = line_construct_pp(&lseg->p[0], &lseg->p[1]);
  1522.     if (p = line_interpt(tmp, line))
  1523.         if (! on_ps(p, lseg)) {
  1524.             PFREE(p);
  1525.             p = NULL;
  1526.         }
  1527.  
  1528.     PFREE(tmp);
  1529.     return(p);
  1530. }
  1531.  
  1532.  
  1533. /*---------------------------------------------------------------------
  1534.  *    close_
  1535.  *        Point of closest proximity between objects.
  1536.  *-------------------------------------------------------------------*/
  1537.  
  1538. /* close_pl - 
  1539.  *    The intersection point of a perpendicular of the line 
  1540.  *    through the point.
  1541.  */
  1542. POINT *
  1543. close_pl(pt, line)
  1544.     POINT    *pt;
  1545.     LINE    *line;
  1546. {
  1547.     POINT    *result;
  1548.     LINE    *tmp;
  1549.     double    invm;
  1550.  
  1551.     result = PALLOCTYPE(POINT);
  1552.     if (FPeq(line->A, -1.0) && FPzero(line->B)) {    /* vertical */
  1553.         result->x = line->C;
  1554.         result->y = pt->y;
  1555.         return(result);
  1556.     } else if (FPzero(line->m)) {            /* horizontal */
  1557.         result->x = pt->x;
  1558.         result->y = line->C;
  1559.         return(result);
  1560.     }
  1561.     /* drop a perpendicular and find the intersection point */
  1562.     invm = -1.0 / line->m;
  1563.     tmp = line_construct_pm(pt, invm);
  1564.     result = line_interpt(tmp, line);
  1565.     return(result);
  1566. }
  1567.  
  1568.  
  1569. /* close_ps - 
  1570.  *    Take the closest endpoint if the point is left, right, 
  1571.  *    above, or below the segment, otherwise find the intersection
  1572.  *    point of the segment and its perpendicular through the point.
  1573.  */
  1574. POINT *
  1575. close_ps(pt, lseg)
  1576.     POINT    *pt;
  1577.     LSEG    *lseg;
  1578. {
  1579.     POINT    *result;
  1580.     LINE    *tmp;
  1581.     double    invm;
  1582.     int    xh, yh;
  1583.  
  1584.     result = NULL;
  1585.     xh = lseg->p[0].x < lseg->p[1].x;
  1586.     yh = lseg->p[0].y < lseg->p[1].y;
  1587.     if (pt->x < lseg->p[!xh].x)
  1588.         result = point_copy(&lseg->p[!xh]);
  1589.     else if (pt->x > lseg->p[xh].x)
  1590.         result = point_copy(&lseg->p[xh]);
  1591.     else if (pt->y < lseg->p[!yh].y)
  1592.         result = point_copy(&lseg->p[!yh]);
  1593.     else if (pt->y > lseg->p[yh].y)
  1594.         result = point_copy(&lseg->p[yh]);
  1595.     if (result)
  1596.         return(result);
  1597.     if (FPeq(lseg->p[0].x, lseg->p[1].x)) {    /* vertical */
  1598.         result->x = lseg->p[0].x;
  1599.         result->y = pt->y;
  1600.         return(result);
  1601.     } else if (FPzero(lseg->m)) {            /* horizontal */
  1602.         result->x = pt->x;
  1603.         result->y = lseg->p[0].y;
  1604.         return(result);
  1605.     }
  1606.  
  1607.     invm = -1.0 / lseg->m;
  1608.     tmp = line_construct_pm(pt, invm);
  1609.     result = interpt_sl(lseg, tmp);
  1610.     return(result);
  1611. }
  1612.  
  1613. /*ARGSUSED*/
  1614. POINT *
  1615. close_pb(pt, box)
  1616.     POINT    *pt;
  1617.     BOX    *box;
  1618. {
  1619.     /* think about this one for a while */
  1620.  
  1621.     return(NULL);
  1622. }
  1623.  
  1624. POINT *
  1625. close_sl(lseg, line)
  1626.     LSEG    *lseg;
  1627.     LINE    *line;
  1628. {
  1629.     POINT    *result;
  1630.     double    *d1, *d2;
  1631.  
  1632.     if (result = interpt_sl(lseg, line))
  1633.         return(result);
  1634.     d1 = dist_pl(&lseg->p[0], line);
  1635.     d2 = dist_pl(&lseg->p[1], line);
  1636.     if (d1 < d2)
  1637.         result = point_copy(&lseg->p[0]);
  1638.     else
  1639.         result = point_copy(&lseg->p[1]);
  1640.  
  1641.     PFREE(d1);
  1642.     PFREE(d2);
  1643.     return(result);
  1644. }
  1645.  
  1646. /*ARGSUSED*/
  1647. POINT *
  1648. close_sb(lseg, box)
  1649.     LSEG    *lseg;
  1650.     BOX    *box;
  1651. {
  1652.     /* think about this one for a while */
  1653.  
  1654.     return(NULL);
  1655. }
  1656.  
  1657. /*ARGSUSED*/
  1658. POINT *
  1659. close_lb(line, box)
  1660.     LINE    *line;
  1661.     BOX    *box;
  1662. {
  1663.     /* think about this one for a while */
  1664.  
  1665.     return(NULL);
  1666. }
  1667.  
  1668. /*---------------------------------------------------------------------
  1669.  *    on_
  1670.  *        Whether one object lies completely within another.
  1671.  *-------------------------------------------------------------------*/
  1672.  
  1673. /* on_pl -
  1674.  *    Does the point satisfy the equation? 
  1675.  */
  1676. long
  1677. on_pl(pt, line)
  1678.     POINT    *pt;
  1679.     LINE    *line;
  1680. {
  1681.     return( FPzero(line->A * pt->x + line->B * pt->y + line->C) );
  1682. }
  1683.  
  1684.  
  1685. /* on_ps -
  1686.  *    Determine colinearity by detecting a triangle inequality.
  1687.  */
  1688. long
  1689. on_ps(pt, lseg)
  1690.     POINT    *pt;
  1691.     LSEG    *lseg;
  1692. {
  1693.     return( point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1])
  1694.            == point_dt(&lseg->p[0], &lseg->p[1]) );
  1695. }
  1696.  
  1697. long
  1698. on_pb(pt, box)
  1699.     POINT    *pt;
  1700.     BOX    *box;
  1701. {
  1702.     return( pt->x <= box->xh && pt->x >= box->xl &&
  1703.            pt->y <= box->yh && pt->y >= box->yl );
  1704. }
  1705.  
  1706. /* on_ppath - 
  1707.  *    Whether a point lies within (on) a polyline.
  1708.  *    If open, we have to (groan) check each segment.
  1709.  *    If closed, we use the old O(n) ray method for point-in-polygon.
  1710.  *        The ray is horizontal, from pt out to the right.
  1711.  *        Each segment that crosses the ray counts as an 
  1712.  *        intersection; note that an endpoint or edge may touch 
  1713.  *        but not cross.
  1714.  *        (we can do p-in-p in lg(n), but it takes preprocessing)
  1715.  */
  1716. #define NEXT(A)    ((A+1) % path->npts)    /* cyclic "i+1" */
  1717.  
  1718. long
  1719. on_ppath(pt, path)
  1720.     POINT    *pt;
  1721.     PATH    *path;
  1722. {
  1723.     int    above, next,    /* is the seg above the ray? */
  1724.     inter,        /* # of times path crosses ray */
  1725.     hi,        /* index inc of higher seg (0,1) */
  1726.     i, n;
  1727.     double a, b, x, 
  1728.     yh, yl, xh, xl;
  1729.  
  1730.     if (! path->closed) {        /*-- OPEN --*/
  1731.         n = path->npts - 1;
  1732.         a = point_dt(pt, &path->p[0]);
  1733.         for (i = 0; i < n; i++) {
  1734.             b = point_dt(pt, &path->p[i+1]);
  1735.             if (FPeq(a+b,
  1736.                      point_dt(&path->p[i], &path->p[i+1])))
  1737.                 return(1);
  1738.             a = b;
  1739.         }
  1740.         return(0);
  1741.     }
  1742.  
  1743.     inter = 0;            /*-- CLOSED --*/
  1744.     above = FPgt(path->p[0].y, pt->y) ? ABOVE : 
  1745.         FPlt(path->p[0].y, pt->y) ? BELOW : UNDEF;
  1746.  
  1747.     for (i = 0; i < path->npts; i++) {
  1748.         hi = path->p[i].y < path->p[NEXT(i)].y;
  1749.         /* must take care of wrap around to original vertex for closed paths */
  1750.         yh = (i+hi < path->npts) ? path->p[i+hi].y : path->p[0].y;
  1751.         yl = (i+!hi < path->npts) ? path->p[i+!hi].y : path->p[0].y;
  1752.         hi = path->p[i].x < path->p[NEXT(i)].x;
  1753.         xh = (i+hi < path->npts) ? path->p[i+hi].x : path->p[0].x;
  1754.         xl = (i+!hi < path->npts) ? path->p[i+!hi].x : path->p[0].x;
  1755.         /* skip seg if it doesn't touch the ray */
  1756.  
  1757.         if (FPeq(yh, yl))    /* horizontal seg? */
  1758.             if (FPge(pt->x, xl) && FPle(pt->x, xh) &&
  1759.                 FPeq(pt->y, yh))
  1760.                 return(1);    /* pt lies on seg */
  1761.             else
  1762.                 continue;    /* skip other hz segs */
  1763.         if (FPlt(yh, pt->y) ||    /* pt is strictly below seg */
  1764.             FPgt(yl, pt->y))    /* strictly above */
  1765.             continue;
  1766.  
  1767.         /* seg touches the ray, find out where */
  1768.  
  1769.         x = FPeq(xh, xl)    /* vertical seg? */
  1770.             ? path->p[i].x    
  1771.                 : (pt->y - path->p[i].y) / 
  1772.                     point_sl(&path->p[i],
  1773.                          &path->p[NEXT(i)]) +
  1774.                              path->p[i].x;
  1775.         if (FPeq(x, pt->x))    /* pt lies on this seg */
  1776.             return(1);
  1777.  
  1778.         /* does the seg actually cross the ray? */
  1779.  
  1780.         next = FPgt(path->p[NEXT(i)].y, pt->y) ? ABOVE : 
  1781.             FPlt(path->p[NEXT(i)].y, pt->y) ? BELOW : above;
  1782.         inter += FPge(x, pt->x) && next != above;
  1783.         above = next;
  1784.     }
  1785.     return(    above == UNDEF ||     /* path is horizontal */
  1786.            inter % 2);        /* odd # of intersections */
  1787. }
  1788.  
  1789.  
  1790. long
  1791. on_sl(lseg, line)
  1792.     LSEG    *lseg;
  1793.     LINE    *line;
  1794. {
  1795.     return( on_pl(&lseg->p[0], line) && on_pl(&lseg->p[1], line) );
  1796. }
  1797.  
  1798. long
  1799. on_sb(lseg, box)
  1800.     LSEG    *lseg;
  1801.     BOX    *box;
  1802. {
  1803.     return( on_pb(&lseg->p[0], box) && on_pb(&lseg->p[1], box) );
  1804. }
  1805.  
  1806. /*---------------------------------------------------------------------
  1807.  *    inter_
  1808.  *        Whether one object intersects another.
  1809.  *-------------------------------------------------------------------*/
  1810.  
  1811. long
  1812. inter_sl(lseg, line)
  1813.     LSEG    *lseg;
  1814.     LINE    *line;
  1815. {
  1816.     POINT    *tmp;
  1817.  
  1818.     if (tmp = interpt_sl(lseg, line)) {
  1819.         PFREE(tmp);
  1820.         return(1);
  1821.     }
  1822.     return(0);
  1823. }
  1824.  
  1825. /*ARGSUSED*/
  1826. long
  1827. inter_sb(lseg, box)
  1828.     LSEG    *lseg;
  1829.     BOX    *box;
  1830. {
  1831.     return(0);
  1832. }
  1833.  
  1834. /*ARGSUSED*/
  1835. long
  1836. inter_lb(line, box)
  1837.     LINE    *line;
  1838.     BOX    *box;
  1839. {
  1840.     return(0);
  1841. }
  1842.  
  1843. /*------------------------------------------------------------------
  1844.  * The following routines define a data type and operator class for
  1845.  * POLYGONS .... Part of which (the polygon's bounding box is built on 
  1846.  * top of the BOX data type.
  1847.  *
  1848.  * make_bound_box - create the bounding box for the input polygon
  1849.  *------------------------------------------------------------------*/
  1850.  
  1851.     /* Maximum number of output digits printed */
  1852. #define P_MAXDIG 12
  1853.  
  1854. /*---------------------------------------------------------------------
  1855.  * Make the smallest bounding box for the given polygon.
  1856.  *---------------------------------------------------------------------*/
  1857. void
  1858. make_bound_box(poly)
  1859. POLYGON *poly;
  1860. {
  1861.     double x1,y1,x2,y2;
  1862.     int npts;
  1863.  
  1864.     npts = poly->npts;
  1865.     x1 = poly_min((double *)poly->pts, npts);
  1866.     x2 = poly_max((double *)poly->pts, npts);
  1867.     y1 = poly_min(((double *)poly->pts)+npts, npts),
  1868.     y2 = poly_max(((double *)poly->pts)+npts, npts);
  1869.     box_fill(&(poly->boundbox), x1, x2, y1, y2); 
  1870. }
  1871.     
  1872. /*------------------------------------------------------------------
  1873.  * polygon_in - read in the polygon from a string specification
  1874.  *              the string is of the form "(f8,f8,f8,f8,...,f8)"
  1875.  *------------------------------------------------------------------*/
  1876. POLYGON *
  1877. poly_in(s)
  1878. char *s;
  1879. {
  1880.     POLYGON *poly;
  1881.     long points;
  1882.     double *xp, *yp, strtod();
  1883.     int i, size;
  1884.  
  1885.     if((points = poly_pt_count(s, ',')) < 0)
  1886.         elog(WARN, "Bad input polyon");
  1887.  
  1888.     size = 2*sizeof(double)*points + sizeof(BOX) + 2*sizeof(int);
  1889.     poly = (POLYGON *)PALLOC(size);
  1890.  
  1891.     if (!PointerIsValid(poly))
  1892.         elog(WARN, "Memory allocation failed, can't input polygon");
  1893.  
  1894.     poly->npts = points;
  1895.     poly->size = size;
  1896.  
  1897.     /* Store all x coords followed by all y coords */
  1898.     xp = (double *) &(poly->pts[0]);
  1899.     yp = (double *) (poly->pts + points*sizeof(double));
  1900.  
  1901.     s++;                /* skip LDELIM */
  1902.  
  1903.     for (i=0; i<points; i++,xp++,yp++)
  1904.     {
  1905.         *xp = strtod(s, &s);
  1906.         s++;                    /* skip delimiter */
  1907.         *yp = strtod(s, &s);
  1908.         s++;                    /* skip delimiter */
  1909.     }
  1910.     make_bound_box(poly);
  1911.     return (poly);
  1912. }
  1913.  
  1914. /*-------------------------------------------------------------
  1915.  * poly_pt_count - count the number of points specified in the
  1916.  *                 polygon.
  1917.  *-------------------------------------------------------------*/
  1918. long
  1919. poly_pt_count(s, delim)
  1920. char *s;
  1921. char delim;
  1922. {
  1923.     long total = 0;
  1924.  
  1925.     if (*s++ != LDELIM)        /* no left delimeter */
  1926.         return (long) -1;
  1927.  
  1928.     while (*s && (*s != RDELIM))
  1929.     {
  1930.         while (*s && (*s != delim))
  1931.             s++;
  1932.         total++;    /* found one */
  1933.         if (*s)
  1934.             s++;    /* bump s past the delimiter */
  1935.     }
  1936.  
  1937.     /* if there was no right delimiter OR an odd number of points */
  1938.  
  1939.     if ((*(s-1) != RDELIM) || ((total%2) != 0))
  1940.         return (long) -1;
  1941.  
  1942.     return (total/2);
  1943. }
  1944.  
  1945. /*---------------------------------------------------------------
  1946.  * poly_out - convert internal POLYGON representation to the 
  1947.  *            character string format "(f8,f8,f8,f8,...f8)"
  1948.  *---------------------------------------------------------------*/
  1949. char *
  1950. poly_out(poly)
  1951. POLYGON *poly;
  1952. {
  1953.     int i;
  1954.     double *xp, *yp;
  1955.     char *output, *outptr;
  1956.  
  1957.     /*-----------------------------------------------------
  1958.      * Get enough space for "(f8,f8,f8,f8,...,f8)"
  1959.      * which P_MAXDIG+1 for each coordinate plus 2
  1960.      * for parens and 1 for the null
  1961.      *-----------------------------------------------------*/
  1962.     output = (char *)PALLOC(2*(P_MAXDIG+1)*poly->npts + 3);
  1963.     outptr = output;
  1964.  
  1965.     if (!output)
  1966.         elog(WARN, "Memory allocation failed, can't output polygon");
  1967.  
  1968.     *outptr++ = LDELIM;
  1969.  
  1970.     xp = (double *) poly->pts;
  1971.     yp = (double *) (poly->pts + (poly->npts * sizeof(double)));
  1972.  
  1973.     sprintf(outptr, "%*g,%*g", P_MAXDIG, *xp++, P_MAXDIG, *yp++);
  1974.     outptr += (2*P_MAXDIG + 1);
  1975.  
  1976.     for (i=1; i<poly->npts; i++,xp++,yp++)
  1977.     {
  1978.         sprintf(outptr, ",%*g,%*g", P_MAXDIG, *xp, P_MAXDIG, *yp);
  1979.         outptr += 2*(P_MAXDIG + 1);
  1980.     }
  1981.     *outptr++ = RDELIM;
  1982.     *outptr = '\0';
  1983.     return (output);
  1984. }
  1985.  
  1986. /*-------------------------------------------------------
  1987.  * Find the largest coordinate out of n coordinates
  1988.  *-------------------------------------------------------*/
  1989. double
  1990. poly_max(coords, ncoords)
  1991. double *coords;
  1992. int ncoords;
  1993. {
  1994.     double max;
  1995.  
  1996.     max = *coords++;
  1997.     ncoords--;
  1998.     while (ncoords--)
  1999.     {
  2000.         if (*coords > max)
  2001.             max = *coords;
  2002.         coords++;
  2003.     }
  2004.     return max;
  2005. }
  2006.  
  2007. /*-------------------------------------------------------
  2008.  * Find the smallest coordinate out of n coordinates
  2009.  *-------------------------------------------------------*/
  2010. double
  2011. poly_min(coords, ncoords)
  2012. double *coords;
  2013. int ncoords;
  2014. {
  2015.     double min;
  2016.  
  2017.     min = *coords++;
  2018.     ncoords--;
  2019.     while (ncoords--)
  2020.     {
  2021.         if (*coords < min)
  2022.             min = *coords;
  2023.         coords++;
  2024.     }
  2025.     return min;
  2026. }
  2027.  
  2028. /*-------------------------------------------------------
  2029.  * Is polygon A strictly left of polygon B? i.e. is
  2030.  * the right most point of A left of the left most point
  2031.  * of B?
  2032.  *-------------------------------------------------------*/
  2033. long
  2034. poly_left(polya, polyb)
  2035. POLYGON *polya, *polyb;
  2036. {
  2037.     double right, left;
  2038.  
  2039.     right = poly_max((double *)polya->pts, polya->npts);
  2040.     left = poly_min((double *)polyb->pts, polyb->npts);
  2041.  
  2042.     return (right < left);
  2043. }
  2044.  
  2045. /*-------------------------------------------------------
  2046.  * Is polygon A overlapping or left of polygon B? i.e. is
  2047.  * the left most point of A left of the right most point
  2048.  * of B?
  2049.  *-------------------------------------------------------*/
  2050. long
  2051. poly_overleft(polya, polyb)
  2052. POLYGON *polya, *polyb;
  2053. {
  2054.     double left, right;
  2055.  
  2056.     left = poly_min((double *)polya->pts, polya->npts);
  2057.     right = poly_max((double *)polyb->pts, polyb->npts);
  2058.  
  2059.     return (left <= right);
  2060. }
  2061.  
  2062. /*-------------------------------------------------------
  2063.  * Is polygon A strictly right of polygon B? i.e. is
  2064.  * the left most point of A right of the right most point
  2065.  * of B?
  2066.  *-------------------------------------------------------*/
  2067. long
  2068. poly_right(polya, polyb)
  2069. POLYGON *polya, *polyb;
  2070. {
  2071.     double right, left;
  2072.  
  2073.     left = poly_min((double *)polya->pts, polya->npts);
  2074.     right = poly_max((double *)polyb->pts, polyb->npts);
  2075.  
  2076.     return (left > right);
  2077. }
  2078.  
  2079. /*-------------------------------------------------------
  2080.  * Is polygon A overlapping or right of polygon B? i.e. is
  2081.  * the right most point of A right of the left most point
  2082.  * of B?
  2083.  *-------------------------------------------------------*/
  2084. long
  2085. poly_overright(polya, polyb)
  2086. POLYGON *polya, *polyb;
  2087. {
  2088.     double right, left;
  2089.  
  2090.     right = poly_max((double *)polya->pts, polya->npts);
  2091.     left = poly_min((double *)polyb->pts, polyb->npts);
  2092.  
  2093.     return (right > left);
  2094. }
  2095.  
  2096. /*-------------------------------------------------------
  2097.  * Is polygon A the same as polygon B? i.e. are all the
  2098.  * points the same?
  2099.  *-------------------------------------------------------*/
  2100. long
  2101. poly_same(polya, polyb)
  2102. POLYGON *polya, *polyb;
  2103. {
  2104.     int i;
  2105.     double *axp, *bxp; /* point to x coordinates for a and b */
  2106.  
  2107.     if (polya->npts != polyb->npts)
  2108.         return 0;
  2109.  
  2110.     axp = (double *)polya->pts;
  2111.     bxp = (double *)polyb->pts;
  2112.  
  2113.     for (i=0; i<polya->npts; axp++, bxp++, i++)
  2114.     {
  2115.         if (*axp != *bxp)
  2116.             return 0;
  2117.     }
  2118.     return 1;
  2119. }
  2120.  
  2121. /*-----------------------------------------------------------------
  2122.  * Determine if polygon A overlaps polygon B by determining if
  2123.  * their bounding boxes overlap.
  2124.  *-----------------------------------------------------------------*/
  2125. long
  2126. poly_overlap(polya, polyb)
  2127. POLYGON *polya, *polyb;
  2128. {
  2129.     return box_overlap(&(polya->boundbox), &(polyb->boundbox));
  2130. }
  2131.  
  2132. /*-----------------------------------------------------------------
  2133.  * Determine if polygon A contains polygon B by determining if A's
  2134.  * bounding box contains B's bounding box.
  2135.  *-----------------------------------------------------------------*/
  2136. long
  2137. poly_contain(polya, polyb)
  2138. POLYGON *polya, *polyb;
  2139. {
  2140.     return box_contain(&(polya->boundbox), &(polyb->boundbox));
  2141. }
  2142.  
  2143. /*-----------------------------------------------------------------
  2144.  * Determine if polygon A is contained by polygon B by determining 
  2145.  * if A's bounding box is contained by B's bounding box.
  2146.  *-----------------------------------------------------------------*/
  2147. long
  2148. poly_contained(polya, polyb)
  2149. POLYGON *polya, *polyb;
  2150. {
  2151.     return box_contained(&(polya->boundbox), &(polyb->boundbox));
  2152. }
  2153.